module Optimisation
   use ObjFunction
   use Equations
   use OptimisationTypes
   use strutils
   use dateutils
   use HTypes
   use Index
   use Bootstrap
   use Retro
   implicit none


contains



!   interface 
!      subroutine OnRoundFinished(iter, delta, best)
!         integer, intent(in) :: iter
!         double precision, intent(in) :: delta
!         double precision, intent(in) :: best
!      end subroutine OnRoundFinished
!   end interface


function RunOpt(Opt)
! Purpose: Sequential search (random sequence) with automated adjustment of increments. 
!          Robust and not so awfully slow.
! In/out:  Delta -
! Out:     jflg  - set to true if the search stops because the max number of iteration 
!                   is reached at min. delta
! Issue:   Test with double precision precision 
! History  When    Who       What
!          ??????  Skagen    Created
!          070521  Skalevik  Translated to F90 from the amci24 optrutine
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   logical :: RunOpt
   double precision :: obj
!   external :: CallBackFunc
   
   integer :: napar
   if(Opt%RunErrorStr /= "") return

   if (associated(Opt%pactive)) then
      napar = size(Opt%pactive)
   else
      napar = 0   
   endif
   
   if ((napar > 0) .and. (Opt%IndexInput%ProcessType /= ptSingleModelRun)) then
      ! Search in the parameter space
      RunOpt = RunGenericSearch(Opt)
   else
      ! single run
      RunOpt = .true.
      obj = objfunc(opt, .true.)
   endif   
   if (RunOpt) then
      RunOpt = OnProcessStageEnd(Opt)
   endif   
end function RunOpt
      
      
function OnProcessStageEnd(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   logical :: OnProcessStageEnd, res
   OnProcessStageEnd = .true.
   ! Process action determination
   select case (Opt%IndexInput%ProcessType)
      case (ptOrdinaryAssessment, ptSingleModelRun)
         call OnOptimisationFinished(Opt)
         return ! No process stage action taken
   end select
   ! Process initiation action
   select case(Opt%Loc%Stage)
      case(0)
         ! Bootstrap process needs inital preparation of data source
         select case (Opt%IndexInput%ProcessType)
            case (ptBootstrap)
               call OnInitiateBootstrap(Opt)
            case (ptRetrospectiveRun)
               call OnInitiateRetro(Opt)
         end select
   end select
   ! Process stage finalization
   call CalculateStockSummary(Opt)
   select case (Opt%IndexInput%ProcessType)
      case (ptBootstrap)
         call OnBootstrapStop(Opt)
      case (ptRetrospectiveRun)
         call OnRetroStop(Opt)
   end select
   if (Opt%Loc%Stage < Opt%IndexInput%NumStages) then
      Opt%Loc%Stage = Opt%Loc%Stage + 1
      ! Process stage reintialization
      select case (Opt%IndexInput%ProcessType)
         case (ptBootstrap)
            call OnBootstrapStart(Opt)
         case (ptRetrospectiveRun)
            call OnRetroStart(Opt)
      end select
      ! After the process stage has been restarted, initialize the run:
      call InitializeRun(Opt)
      OnProcessStageEnd = .false.

      !call SavePopulationModel(Opt)
      !call system("pause")

   else
      ! Process is finished -> reload run again 
      res = LoadRun(Opt)
   endif
end function OnProcessStageEnd
 
 
subroutine WriteOutAllParameters(Opt, filename)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   character(len=*) :: filename
   integer :: NFleets, NYears, NAges, Cat, iFleet, iYear, iAge
   type(TParamElm), pointer :: pParamElm
   type(TParamCube), pointer :: pParamCube
   
   open(159, file = filename, status = 'replace')
   do Cat = cFirstCategory, cLastCategory
      pParamCube => Opt%ParamCubes(Cat)
      NFleets = size(pParamCube%Elm, 1)
      do iFleet = 0, NFleets - 1
         NYears = size(pParamCube%Elm, 2)
         do iYear = 0, NYears - 1
            NAges = size(pParamCube%Elm, 3)
            do iAge = 0, NAges - 1
               pParamElm => GetParamElm(Opt, Cat, iFleet, iYear, iAge, .false.)
               write(159,*) Cat, iFleet, iYear, iAge, pParamElm%Value
            enddo
         enddo
      enddo
   enddo
   close(159)
end subroutine WriteOutAllParameters


subroutine ImproveParameters(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: napar, i, j, j1, jbest, jpar, Ip, IPar, k
   double precision :: pvalue, s1, obj, objbefore 
   double precision :: seqj
   logical :: TestError = .false.
!   double precision, pointer  :: activebefore(:)
 !  integer :: iab
   logical :: ErrorDetected
   
   Opt%Loc%Improved = .false.
   ErrorDetected = .false.
   !if (TestError) then
   !   allocate(activebefore(0 : size(Opt%pactive) - 1))
   !endif
   ! Order parameters random with fixed seed
   napar = size(Opt%pactive)
   do i = 1, napar
      seqj = rndn(Opt%seed)
      Opt%seq(i) = seqj
   end do
   do i = 1, napar
      s1 = 1.1
      j1 = i
      do j = 1, napar
         seqj = Opt%seq(j)
         if ((seqj < s1) .and. (0 < seqj)) then
            s1 = seqj
            j1 = j
         end if
      end do
      Opt%iseq(i) = j1
      Opt%seq(j1) = 0.0
   end do
   ! Try each parameter up or down and test if improved
   if(Opt%DumpTuningDetails) then
      open(cTuningFile, file = GetWorkFileName(Opt, s_TuningFile), action = 'write', position =  'append')
   endif
   do Ip = 1, napar 
      Ipar = Opt%iseq(Ip) - 1
      ! pactive is the list of the active parameter value
      pvalue = Opt%pactive(Ipar)%Value
      jbest = 0
      objbefore = Opt%Loc%ObjFunc
      
      !if (TestError) then
      !   do iab = 0, size(activebefore) - 1
      !      activebefore(iab) = Opt%pactive(iab)%value
      !   enddo
      !endif
      !call WriteOutAllParameters(Opt, "c:\before.txt")
      ! jbest holds the best alternative; 0 unless there is improvement
      do jpar = -1, 1 
         ! Aasmund 12.06.07 - Added step=2 to avoid objfunc() in the case jpar=0 (current obj)
         call SetPValueEx(Opt%pactive(Ipar), (1 + Opt%Loc%Delta * jpar) * pvalue)
         obj = objfunc(Opt, .false.)
         
         if ((obj < Opt%Loc%ObjFunc) .and. (abs(obj - Opt%Loc%ObjFunc) > 0.000001)) then
            jbest = jpar
            Opt%Loc%ObjFunc = obj
            ! Improved signals that the value of this parameter was improved
            Opt%Loc%Improved = .true.
         end if
      end do
      ! Ensure that pvector holds the right value when leaving this parameter
      pvalue = (1 + Opt%Loc%Delta * jbest) * pvalue
      call SetPValueEx(Opt%pactive(Ipar), pvalue)
      
      if (TestError .and. (jbest == 0)) then
         ! Check everything
         obj = objfunc(Opt, .false.)
         if (obj /= objbefore) then
            !write(*,*) "Error after"
            !call WriteOutAllParameters(Opt, "c:\after.txt")
            !call system("pause")
            !write(*,*) "ERROR Detected"
            !do iab = 0, size(activebefore) - 1
!           !    write (*,* ) activebefore(iab), Opt%pactive(iab)%value
            !   if (activebefore(iab) /= Opt%pactive(iab)%value) then
            !      call MessageOut(Opt, "Error: Parameter " // lrtrim(inttostr(iab)) // "corrupted by objfunc " // &
            !         " Cat:" // lrtrim(inttostr(Opt%pactive(iab)%Cat)) // & 
            !         " Index:" // lrtrim(inttostr(Opt%pactive(iab)%Index)) // & 
            !         " Year:" // lrtrim(inttostr(Opt%pactive(iab)%Index)) // & 
            !         " Age:" // lrtrim(inttostr(Opt%pactive(iab)%Index))) 
            !   endif
            !enddo
         endif
      endif
      if(Opt%DumpTuningDetails .and. (jbest /= 0)) then
         write(cTuningFile, *) IPar, Opt%Loc%Delta * jbest, pvalue
      endif
   end do
   if(Opt%DumpTuningDetails) then
      close(cTuningFile)
   endif
!   if (TestError) then
!      deallocate(activebefore)
!   endif
end subroutine ImproveParameters


subroutine OnOptimisationFinished(Opt)
! Purpose: Do what is necessary when the optimization is finished
!          1. Calculate the object function. Reason:
!              - produce unweighted residuals around the model fit.ensure the rescaled residuals. 
!          2. Normalize(scale) the selection and year factor to the reference age range
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   ! Calculate residuals and stream out phi values (Finished parameter = true)
   ! Save population model (N and F), summary, datasource residuals, and active parameters
   call SavePopulationModel(Opt)
   ! Save tuned model parameters to out file.
   call SaveParameters(Opt, .true.)
end subroutine OnOptimisationFinished


function RunGenericSearch(Opt)
! Purpose: Sequential search (random sequence) with automated adjustment of increments. 
!          Robust and not so awfully slow.
! In/out:  Delta -
! Out:     jflg  - set to true if the search stops because the max number of iteration 
!                   is reached at min. delta
! Issue:   Test with double precision precision 
! History  When    Who       What
!          ??????  Skagen    Created
!          070521  Skalevik  Translated to F90 from the amci24 optrutine

   implicit none
   type(TOptimisation), intent(inout) :: Opt
   logical :: RunGenericSearch
   logical :: ChangeDelta
   !double precision, pointer :: pobj(:)
   integer, parameter :: c_MaxIterations = 25
   double precision, parameter :: c_MinDelta = 0.0001d0
   double precision, parameter :: c_MaxDelta = 0.4d0
   double precision, parameter :: c_DeltaIncreaseRate = 2.0d0
   double precision, parameter :: c_DeltaDecreaseRate = 0.4d0
   double precision :: obj
   
   !allocate(pobj(size(Opt%pactive)))
   ! Start with increment delta (=0.25), then delta = delta / 5
   ! c_MaxIterations is the max number of rounds taken before
   ! an increase in delta is made. The choice has some impact
   ! on the convergence rate
   ! Actions to do 
   ! - Prepare delta and iteration for this round
   ! - Determine if finished based on Improved state  and iteration
   ! - Improve parameters and update improved state.
   ! Change delta on iteration 25 or if no improvement detected
   obj = objfunc(Opt, .false.)
   ChangeDelta = (Opt%Loc%Iter == c_MaxIterations) .or. (.not. Opt%Loc%Improved)
   if (ChangeDelta) then
      if (c_MinDelta < Opt%Loc%Delta) then
         ! If last iteration is processed, double delta
         if (Opt%Loc%Iter == c_MaxIterations) then
            if (Opt%Loc%Delta < c_MaxDelta) Opt%Loc%Delta = Opt%Loc%Delta * c_DeltaIncreaseRate
         else 
            Opt%Loc%Delta = Opt%Loc%Delta * c_DeltaDecreaseRate
         end if
         !if (Opt%DebugMode) then
            !call OnOptimisationFinished(Opt)
            call SaveParameters(Opt, .true.)
         !end if
      else
         RunGenericSearch = .true.
         return
      end if
      Opt%Loc%Iter = 1
   else
      Opt%Loc%Iter = Opt%Loc%Iter + 1
   end if

   ! Improvement check:
   call ImproveParameters(Opt)
   obj = objfunc(Opt, .false.)
   RunGenericSearch = .false.
end function RunGenericSearch


end module Optimisation     